#############################################
# differential abundance analysis
#############################################

library(ggordiplots)
library(microViz)
library("vegan")
library(phyloseqCompanion)
library(ggplot2)
library(Rmisc)
library(ggsignif)
library(ggpubr)
library(dplyr)
library("gridExtra")
library(microbiomeMarker)
library(tidyr)

# set scientific notation off
options(scipen=999)

# set directory
setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code")

# load phyloseq variables
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/phyloseq final variables.RData")

# need to normalize in some way before differential analysis
# can do rarefy, use relative abundance

# transform using log transformation and pseudocount of 1 (don't transform unless otherwise indicated)
#genus.tr <- tax_transform(genus.glom, trans = "log", rank = "Genus", zero_replace = 1)

# read in soil data
soil.data <- read.csv("soil.data.final.csv")

# add covariates to phyloseq object
sample_data(physeq)$organic.soil.C.N <- soil.data$organic.soil.C.N[1:96]
sample_data(physeq)$mineral.soil.C.N <- soil.data$mineral.soil.C.N[1:96]
sample_data(physeq)$organic.soil.GWC <- soil.data$organic.soil.GWC[1:96]
sample_data(physeq)$carcass.deposition <- soil.data$carcass.deposition[1:96]
sample_data(physeq)$ph <- soil.data$ph[1:96]
sample_data(physeq)$decomposing.carcass <- soil.data$decomposing.carcass[1:96]
sample_data(physeq)$white.spruce <- soil.data$white.spruce[1:96]
sample_data(physeq)$paper.birch <- soil.data$paper.birch[1:96]
sample_data(physeq)$dwarf.birch <- soil.data$dwarf.birch[1:96]
sample_data(physeq)$moss <- soil.data$moss[1:96]
sample_data(physeq)$vaccinium <- soil.data$vaccinium[1:96]
sample_data(physeq)$fern <- soil.data$fern[1:96]
sample_data(physeq)$kinnikinnick <- soil.data$kinnikinnick
sample_data(physeq)$heather <- soil.data$heather[1:96]
sample_data(physeq)$eriophorum <- soil.data$eriophorum[1:96]
sample_data(physeq)$fireweed <- soil.data$fireweed[1:96]
sample_data(physeq)$horsetail <- soil.data$horsetail[1:96]
sample_data(physeq)$alder <- soil.data$alder[1:96]
sample_data(physeq)$hogsweed <- soil.data$hogsweed[1:96]
sample_data(physeq)$willow <- soil.data$willow[1:96]
sample_data(physeq)$redberry <- soil.data$redberry[1:96]
sample_data(physeq)$azalea <- soil.data$azalea[1:96]
sample_data(physeq)$slope <- soil.data$slope[1:96]

# use relative abundance
physeq.prop <- transform_sample_counts(physeq, function(x) x/sum(x))

#################################################
# Wilcoxon test to compare the abundance of different genera by Hansen close bank
###################################################

# subset by Hansen and agglomerate by genus
physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
physeq_hansen_close <- subset_samples(physeq_hansen, Distance < 7)
genus.glom <- tax_glom(physeq.prop, taxrank = 'Genus')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Genus <- as.character(dat_genus$Genus)
dat_genus$Genus <- dat_genus$Genus %>% replace_na('Unclassified')
genus.abund <- subset(dat_genus, select = c("Sample", "Abundance", "Stream", "Transect", "Side", "Distance", "organic.soil.C.N", 
    "mineral.soil.C.N", "organic.soil.GWC","carcass.deposition", "ph", "decomposing.carcass", "Genus"))

# create unique number for each genera
genus.abund <- genus.abund %>%
  group_by(Genus) %>%
  mutate(ID = cur_group_id())

# create unique number for each sample
genus.abund <- genus.abund %>%
  group_by(Sample) %>%
  mutate(sample.ID = cur_group_id())

genus.abund <- as.data.frame(genus.abund)

# remove S4 SR 1, 3, and 6 (if not using decomposing carcass variable)
genus.abund <- genus.abund[-which(genus.abund$Stream=="Hansen" & genus.abund$Distance=="1" & genus.abund$Side=="SR" & genus.abund$Transect=="S4"),]
genus.abund <- genus.abund[-which(genus.abund$Stream=="Hansen" & genus.abund$Distance=="3" & genus.abund$Side=="SR" & genus.abund$Transect=="S4"),]
genus.abund <- genus.abund[-which(genus.abund$Stream=="Hansen" & genus.abund$Distance=="6" & genus.abund$Side=="SR" & genus.abund$Transect=="S4"),]

# create vector for each genera and p value
pvalue <- c(1:308)

for (i in 1:308){
  genus.final <- genus.abund[which(genus.abund$ID==i),]
  left <- genus.final$Abundance[which(genus.final$Side=="SL")]
  right <- genus.final$Abundance[which(genus.final$Side=="SR")]
  w <- wilcox.test(left, right, exact=F)
  pvalue[i] <- w$p.value
}

# create dataframe with genus and p value combined
genus.p <- data.frame(levels(factor(genus.abund$Genus)), pvalue)
genus.p <- genus.p[order(pvalue),]
genus.sig <- genus.p[which(genus.p$pvalue < 0.051),]
colnames(genus.sig) <- c("Genus","P value")
View(genus.sig)
write.csv(genus.sig, "Sig genus abun by Hansen 1-6 side.csv")

# export means by stream for each genus to make bar graph in excel
final <- genus.abund[genus.abund$Genus %in% genus.sig$Genus,]
tgc <- summarySE(final, measurevar="Abundance", groupvars=c("Genus"))
tgc_order <- tgc[order(-tgc$Abundance),]
tgc_order
#write.csv(tgc, "means of sig genus abundance by Hansen side.csv")

# do each genus manually
# choose genus
genus.final <- genus.abund[which(genus.abund$Genus=="Paxillus"),]
left <- genus.final$Abundance[which(genus.final$Side=="SL")]
right <- genus.final$Abundance[which(genus.final$Side=="SR")]

# then use pairwise comparisons between group levels with corrections for multiple testing
# Shapiro-Wilk normality test for the differences (if p > 0.05, can assume normality)
d <- with(genus.final, 
          Abundance[Side == "SL"] - Abundance[Side == "SR"])
shapiro.test(d) 

# use either parametric t-test (assumes normality) or non-parametric Wilcoxon test
t.test(left,right, paired=TRUE)
wilcox.test(left, right, paired=TRUE, exact=F)

# plot a fungal genera
fungus <- genus.abund[which(genus.abund$Genus=="Paxillus"),]
tgc <- summarySE(fungus, measurevar="Abundance", groupvars=c("Side"))
ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", 
   stat = "summary", fun = "mean") + theme_classic()+ 
  theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
                          axis.text.y = element_text(size = 14), legend.position="none") +
  ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

Paxillus <- genus.abund[which(genus.abund$Genus=="Paxillus"),]
Boletus <- genus.abund[which(genus.abund$Genus=="Boletus"),]
Cortinarius <- genus.abund[which(genus.abund$Genus=="Cortinarius"),]
Pseudotomentella <- genus.abund[which(genus.abund$Genus=="Pseudotomentella"),]
Inocybe <- genus.abund[which(genus.abund$Genus=="Inocybe"),]
Clavulina <- genus.abund[which(genus.abund$Genus=="Clavulina"),]
Mortierella <- genus.abund[which(genus.abund$Genus=="Mortierella"),]
Oidiodendron <- genus.abund[which(genus.abund$Genus=="Oidiodendron"),]

tgc <- summarySE(Paxillus, measurevar="Abundance", groupvars=c("Side"))
paxillus <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Paxillus")+
  theme_classic() + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.title.y = element_text(size = 14), 
                          axis.text.y = element_text(size = 14), legend.position="none") + ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9))

tgc <- summarySE(Boletus, measurevar="Abundance", groupvars=c("Side"))
tgc[1,5] <- 0.00001
boletus <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", 
  stat = "summary", fun = "mean") + theme_classic() + ggtitle("Boletus")+
  theme_classic() + theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.title.y = element_blank(), 
                          axis.text.y = element_text(size = 14), legend.position="none") + ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9))
boletus

tgc <- summarySE(Cortinarius, measurevar="Abundance", groupvars=c("Side"))
cortinarius <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Cortinarius")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none")

tgc <- summarySE(Pseudotomentella, measurevar="Abundance", groupvars=c("Side"))
pseudotomentella <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Pseudoomentella")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none")

tgc <- summarySE(Inocybe, measurevar="Abundance", groupvars=c("Side"))
inocybe <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Inocybe")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none")

tgc <- summarySE(Clavulina, measurevar="Abundance", groupvars=c("Side"))
clavulina <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Clavulina")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none")

tgc <- summarySE(Mortierella, measurevar="Abundance", groupvars=c("Side"))
mortierella <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Mortierella")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.ticks.x=element_blank(), axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_x_discrete(labels=c("Salmon \n enhanced", "Salmon \n depleted"))

tgc <- summarySE(Oidiodendron, measurevar="Abundance", groupvars=c("Side"))
oidiodendron <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Oidiodendron")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.ticks.x=element_blank(), axis.text.y = element_text(size = 14),axis.title.y = element_blank(),
        legend.position="none") + ylab("Relative Abundance") + scale_x_discrete(labels=c("Salmon \n enhanced", "Salmon \n depleted")) + 
  geom_signif(comparisons = list(c("SL", "SR")), 
              map_signif_level=TRUE)

p <- grid.arrange(paxillus, boletus, cortinarius, pseudotomentella, inocybe, clavulina, mortierella, oidiodendron, nrow=4, ncol=2)
p

setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing")
ggsave(plot=p, width=7, height=10, dpi=300, filename = "Fig4.jpg")

#################################################
# Wilcoxon test to compare the abundance of different families by Hansen close bank
###################################################

# subset by Hansen and agglomerate by genus
physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
physeq_hansen_close <- subset_samples(physeq_hansen, Distance < 7)
genus.glom <- tax_glom(physeq_hansen_close, taxrank = 'Family')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Family <- as.character(dat_genus$Family)
dat_genus$Family <- dat_genus$Family %>% replace_na('Unclassified')
genus.abund <- dat_genus[,c(3,5,7,13)]

# create unique number for each genera
genus.abund <- genus.abund %>%
  group_by(Family) %>%
  mutate(ID = cur_group_id())

genus.abund <- as.data.frame(genus.abund)

# create vector for each genera and p value
pvalue <- c(1:154)

for (i in 1:154){
  genus.final <- genus.abund[which(genus.abund$ID==i),]
  left <- genus.final$Abundance[which(genus.final$Side=="SL")]
  right <- genus.final$Abundance[which(genus.final$Side=="SR")]
  w <- wilcox.test(left, right, exact=F)
  pvalue[i] <- w$p.value
}

# create dataframe with genus and p value combined
genus.p <- data.frame(levels(factor(genus.abund$Family)), pvalue)
genus.p <- genus.p[order(pvalue),]
genus.sig <- genus.p[which(genus.p$pvalue < 0.051),]
colnames(genus.sig) <- c("Genus","P value")
View(genus.sig)
write.csv(genus.sig, "Sig family abun by Hansen 1-6 side.csv")

# choose genus
genus.final <- genus.abund[which(genus.abund$Family=="Russulaceae"),]
left <- genus.final$Abundance[which(genus.final$Side=="SL")]
right <- genus.final$Abundance[which(genus.final$Side=="SR")]

# there are four families with significant differences between sides

Inocybaceae <- genus.abund[which(genus.abund$Family=="Inocybaceae"),]
Ophiocordycipitaceae <- genus.abund[which(genus.abund$Family=="Ophiocordycipitaceae"),]
Russulaceae <- genus.abund[which(genus.abund$Family=="Russulaceae"),]
Pezizomycotina_fam_Incertae_sedis <- genus.abund[which(genus.abund$Family=="Pezizomycotina_fam_Incertae_sedis"),]
Myxotrichaceae <- genus.abund[which(genus.abund$Family=="Myxotrichaceae"),]
Cortinariaceae <- genus.abund[which(genus.abund$Family=="Cortinariaceae"),]

tgc <- summarySE(Inocybaceae, measurevar="Abundance", groupvars=c("Side"))
Inocybaceae <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + 
  ggtitle("Inocybaceae")+ theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
  axis.text.y = element_text(size = 14), legend.position="none", plot.title = element_text(size = 10)) + ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) 

tgc <- summarySE(Ophiocordycipitaceae, measurevar="Abundance", groupvars=c("Side"))
Ophiocordycipitaceae <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Ophiocordycipitaceae")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none", plot.title = element_text(size = 10))

tgc <- summarySE(Russulaceae, measurevar="Abundance", groupvars=c("Side"))
Russulaceae <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Russulaceae")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none", plot.title = element_text(size = 10))

tgc <- summarySE(Pezizomycotina_fam_Incertae_sedis, measurevar="Abundance", groupvars=c("Side"))
Pezizomycotina_fam_Incertae_sedis <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + 
  ggtitle("Pezizomycotina \n Incertae_sedis")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none", plot.title = element_text(size = 10))

tgc <- summarySE(Myxotrichaceae, measurevar="Abundance", groupvars=c("Side"))
Myxotrichaceae <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + 
  ggtitle("Myxotrichaceae")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none", plot.title = element_text(size = 10))

tgc <- summarySE(Cortinariaceae, measurevar="Abundance", groupvars=c("Side"))
Cortinariaceae <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + 
  ggtitle("Cortinariaceae")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_blank(), axis.text.y = element_text(size = 14),
        legend.position="none", plot.title = element_text(size = 10))

p <- grid.arrange(Inocybaceae, Ophiocordycipitaceae, Russulaceae, Pezizomycotina_fam_Incertae_sedis, Cortinariaceae, Myxotrichaceae, nrow=1, ncol=6)
p

setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing")
ggsave(plot=p, width=10, height=3, dpi=300, filename = "Fig1b.jpg")

##################################################################################################
# Wilcoxon test to compare the abundance of different species by Hansen close bank
####################################################################################################

# use relative abundance
physeq.prop <- transform_sample_counts(physeq, function(x) x/sum(x))

# subset samples
physeq_nobog <- subset_samples(physeq.prop, Transect != "Bog")

# subset by Hansen and agglomerate by genus
physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
physeq_hansen_close <- subset_samples(physeq_hansen, Distance < 7)
species.glom <- tax_glom(physeq_nobog, taxrank = 'Species')

# get species abundance matrix from phyloseq object
dat_species <- psmelt(species.glom)
dat_species$Species <- as.character(dat_species$Species)
dat_species$Species <- dat_species$Species %>% replace_na('Unclassified')
species.abund <- subset(dat_species, select = c("Sample", "Abundance", "Stream", "Transect", "Side", "Distance", "organic.soil.C.N", "organic.soil.GWC",
                                                "carcass.deposition", "ph", "decomposing.carcass", "Genus", "Species"))

# create unique number for each species
species.abund <- species.abund %>%
  group_by(Species) %>%
  mutate(ID = cur_group_id())

# create unique number for each sample
species.abund <- species.abund %>%
  group_by(Sample) %>%
  mutate(sample.ID = cur_group_id())

species.abund <- as.data.frame(species.abund)

# remove S4 SR 1, 3, and 6
species.abund <- species.abund[-which(species.abund$Stream=="Hansen" & species.abund$Distance=="1" & species.abund$Side=="SR" & species.abund$Transect=="S4"),]
species.abund <- species.abund[-which(species.abund$Stream=="Hansen" & species.abund$Distance=="3" & species.abund$Side=="SR" & species.abund$Transect=="S4"),]
species.abund <- species.abund[-which(species.abund$Stream=="Hansen" & species.abund$Distance=="6" & species.abund$Side=="SR" & species.abund$Transect=="S4"),]

# create vector for each genera and p value
pvalue <- c(1:531)

for (i in 1:531){
  species.final <- species.abund[which(species.abund$ID==i),]
  left <- species.final$Abundance[which(species.final$Side=="SL")]
  right <- species.final$Abundance[which(species.final$Side=="SR")]
  w <- wilcox.test(left, right, exact=F)
  pvalue[i] <- w$p.value
}

# create dataframe with species and p value combined
species.p <- data.frame(levels(factor(species.abund$Species)), pvalue)
species.p <- species.p[order(pvalue),]
species.sig <- species.p[which(species.p$pvalue < 0.051),]
colnames(species.sig) <- c("species","P value")
View(species.sig)
write.csv(species.sig, "Sig species abun by Hansen 1-6 side.csv")

# export means by stream for each species to make bar graph in excel
final <- species.abund[species.abund$species %in% species.sig$species,]
tgc <- summarySE(final, measurevar="Abundance", groupvars=c("species"))
tgc_order <- tgc[order(-tgc$Abundance),]
tgc_order
#write.csv(tgc, "means of sig species abundance by Hansen side.csv")

# do each species manually
# choose species
species.final <- species.abund[which(species.abund$Species=="involutus"),]
left <- species.final$Abundance[which(species.final$Side=="SL")]
right <- species.final$Abundance[which(species.final$Side=="SR")]

# then use pairwise comparisons between group levels with corrections for multiple testing
# Shapiro-Wilk normality test for the differences (if p > 0.05, can assume normality)
d <- with(species.final, 
          Abundance[Side == "SL"] - Abundance[Side == "SR"])
shapiro.test(d) 

# use either parametric t-test (assumes normality) or non-parametric Wilcoxon test
t.test(left,right)
wilcox.test(left, right, exact=F)

# plot long distance species by close side
fungus <- species.abund[which(species.abund$Species=="involutus"),]
tgc <- summarySE(fungus, measurevar="Abundance", groupvars=c("Side"))
involutus <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic()+
  theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
  axis.text.y = element_text(size = 14),legend.position="none") +ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9))  + ggtitle("Paxillus involutus") + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

fungus <- species.abund[which(species.abund$Species=="corsicus"),]
tgc <- summarySE(fungus, measurevar="Abundance", groupvars=c("Side"))
corsicus <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic()+
  theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
                          axis.text.y = element_text(size = 14),legend.position="none") +ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + ggtitle("Alpova corsicus") + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

fungus <- species.abund[which(species.abund$Species=="edulis"),]
tgc <- summarySE(fungus, measurevar="Abundance", groupvars=c("Side"))
edulis <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic()+
  theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
                          axis.text.y = element_text(size = 14),legend.position="none") +ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9))  + ggtitle("Boletus edulis") + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

p <- grid.arrange(involutus, corsicus, edulis, xerocomus, nrow=1, ncol=4)
p

# plot significantly different species by close side

setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing")
ggsave(plot=p, width=10, height=4, dpi=300, filename = "species by bank.jpg")

napipes <- species.abund[which(species.abund$Species=="napipes"),]
inflatum <- species.abund[which(species.abund$Species=="inflatum"),]
truncatum <- species.abund[which(species.abund$Species=="truncatum"),]
cystojenkinii <- species.abund[which(species.abund$Species=="cystojenkinii"),]
vietus <- species.abund[which(species.abund$Species=="vietus"),]
botrytis <- species.abund[which(species.abund$Species=="botrytis"),]
sublilacina <- species.abund[which(species.abund$Species=="sublilacina"),]

tgc <- summarySE(napipes, measurevar="Abundance", groupvars=c("Side"))
napipes <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", 
           stat = "summary", fun = "mean") + theme_classic() + ggtitle("Inocybe napipes")+
  theme_classic() + theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), 
  axis.text.y = element_text(size = 14), legend.position="none") + ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(inflatum, measurevar="Abundance", groupvars=c("Side"))
inflatum <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Ophiocordycipitaceae \n inflatum")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14), axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(truncatum, measurevar="Abundance", groupvars=c("Side"))
truncatum <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Oidiodendron truncatum")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(cystojenkinii, measurevar="Abundance", groupvars=c("Side"))
cystojenkinii <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Mortierella \n cystojenkinii")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(vietus, measurevar="Abundance", groupvars=c("Side"))
vietus <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Lactarius vietus")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.text.y = element_text(size = 14),axis.title.y = element_text(size = 14),
        legend.position="none") + ylab("Relative Abundance") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(botrytis, measurevar="Abundance", groupvars=c("Side"))
botrytis <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Amblyosporium botrytis")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

tgc <- summarySE(sublilacina, measurevar="Abundance", groupvars=c("Side"))
sublilacina <- ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Tomentella sublilacina")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.x=element_text(size = 14),axis.text.x=element_text(size = 14),axis.title.y = element_text(size = 14), axis.text.y = element_text(size = 14),
        legend.position="none") + scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

p <- grid.arrange(napipes, inflatum, truncatum, cystojenkinii, vietus, botrytis, sublilacina, nrow=2, ncol=4)
p

setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing")
ggsave(plot=p, width=10, height=6, dpi=300, filename = "indicator species by close bank.jpg")

#################################################
# Kruskal Wallis test to compare the abundance of different genera between streams
###################################################


results <- run_lefse(physeq.prop,
                     wilcoxon_cutoff = 0.1,
                     norm = "none",
                     group = "Stream",
                     kw_cutoff = 0.1,
                     multigrp_strat = T,
                     lda_cutoff = 2
)

# agglomerate by genus
genus.glom <- tax_glom(physeq.prop, taxrank = 'Genus')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Genus <- as.character(dat_genus$Genus)
dat_genus$Genus <- dat_genus$Genus %>% replace_na('Unclassified')
genus.abund <- dat_genus[,c(3,5,14)]

# create unique number for each genera
genus.abund <- transform(genus.abund,id=as.numeric(factor(Genus)))

# create vector for each genera and p value
pvalue <- c(1:308)
w1 <- c(1:308)
w2 <- c(1:308)
w3 <- c(1:308)

for (i in 1:308){
  genus.final <- genus.abund[which(genus.abund$id==i),]
  k <- kruskal.test(Abundance ~ Stream, data = genus.final)
  pvalue[i] <- k$p.value
  w <- pairwise.wilcox.test(genus.final$Abundance, genus.final$Stream, 
                            p.adjust.method = "BH", exact=F)
  w1[i] <- w$p.value[1,1]
  w2[i] <- w$p.value[2,1]
  w3[i] <- w$p.value[2,2]
}

# create dataframe with genus and p value combined
genus.p <- data.frame(levels(factor(genus.abund$Genus)), pvalue, w1, w2, w3)
genus.p <- genus.p[order(pvalue),]
genus.sig <- genus.p[which(genus.p$pvalue < 0.055),]
colnames(genus.sig) <- c("Genus","P value","Hansen_Happy","Hansen_Yako","Happy_Yako")
genus.sig <- genus.sig[,-2]
View(genus.sig)
write.csv(genus.sig, "Sig genus abun by stream.csv")

# export means by stream for each genus to make bar graph in excel
final <- genus.abund[genus.abund$Genus %in% genus.sig$Genus,]
tgc <- summarySE(final, measurevar="Abundance", groupvars=c("Genus", "Stream"))
write.csv(tgc, "means of sig genus abundance by stream.csv")

# do each genus manually
# choose genus
genus.final <- genus.abund[which(genus.abund$Genus=="Amanita"),]

# but really need to run Kruskas Wallis test for comparing more than two groups in order
# not to increase Type I error of multiple t-tests
kruskal.test(Abundance ~ Stream, data = genus.final)

# then use pairwise comparisons between group levels with corrections for multiple testing
pairwise.wilcox.test(genus.final$Abundance, genus.final$Stream, p.adjust.method = "BH",
                     exact=F)

tgc <- summarySE(genus.final, measurevar="Abundance", groupvars=c("Stream"))
ggplot(tgc, aes(Stream, Abundance, fill=Stream)) + geom_bar(position = "dodge", 
                                                            stat = "summary", fun = "mean") + theme_classic() + ggtitle("Amanita")+
  theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
         axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
         axis.title= element_text(face = "bold", size = 14, colour = "black"), 
         panel.background = element_blank(), legend.position="none",
         panel.border = element_rect(fill = NA, colour = "black")) + 
  ylab("Relative Abundance") +
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) 

#################################################
# Kruskal Wallis test to compare the abundance of different genera by transect at Hansen
# Happy, and Yako
###################################################

# agglomerate by genus
physeq_nobog <- subset_samples(physeq.prop, Transect != "Bog")
physeq_yako <- subset_samples(physeq.prop, Stream == "Yako")

genus.glom <- tax_glom(physeq_yako, taxrank = 'Genus')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Genus <- as.character(dat_genus$Genus)
dat_genus$Genus <- dat_genus$Genus %>% replace_na('Unclassified')
genus.abund <- dat_genus[,c(3,5,6,14)]

# create unique number for each genera
genus.abund <- transform(genus.abund,id=as.numeric(factor(Genus)))

# create vector for each genera and p value
pvalue <- c(1:308)
w1 <- c(1:308)
w2 <- c(1:308)
w3 <- c(1:308)

for (i in 1:308){
  genus.final <- genus.abund[which(genus.abund$id==i),]
  k <- kruskal.test(Abundance ~ Transect, data = genus.final)
  pvalue[i] <- k$p.value
  w <- pairwise.wilcox.test(genus.final$Abundance, genus.final$Transect, 
                            p.adjust.method = "BH", exact=F)
  w1[i] <- w$p.value[1,1]
  w2[i] <- w$p.value[2,1]
  w3[i] <- w$p.value[2,2]
}

# create dataframe with genus and p value combined
genus.p <- data.frame(levels(factor(genus.abund$Genus)), pvalue, w1, w2, w3)
genus.p <- genus.p[order(pvalue),]
genus.sig <- genus.p[which(genus.p$pvalue < 0.051),]
#colnames(genus.sig) <- c("Genus","P value","S1_S4","S1_S6","S4_S6")
#colnames(genus.sig) <- c("Genus","P value","H1_H2","H1_H3","H2_H3")
colnames(genus.sig) <- c("Genus","P value","Y1_Y2","Y1_Y3","Y2_Y3")
genus.sig <- genus.sig[,-2]
View(genus.sig)
#write.csv(genus.sig, "Sig genus abun by stream.csv")

# export means by stream for each genus to make bar graph in excel
final <- genus.abund[genus.abund$Genus %in% genus.sig$Genus,]
tgc <- summarySE(final, measurevar="Abundance", groupvars=c("Genus"))
tgc_order <- tgc[order(-tgc$Abundance),]
tgc_order

# sort means by abundance in order to only 

# do each genus manually
# choose genus
genus.final <- genus.abund[which(genus.abund$Genus=="Lachnellula"),]

# but really need to run Kruskas Wallis test for comparing more than two groups in order
# not to increase Type I error of multiple t-tests
kruskal.test(Abundance ~ Transect, data = genus.final)

# then use pairwise comparisons between group levels with corrections for multiple testing
pairwise.wilcox.test(genus.final$Abundance, genus.final$Transect, p.adjust.method = "BH",
                     exact=F)

tgc <- summarySE(genus.final, measurevar="Abundance", groupvars=c("Transect"))
ggplot(tgc, aes(Transect, Abundance, fill=Transect)) + geom_bar(position = "dodge", 
  stat = "summary", fun = "mean") + theme_classic() + ggtitle("Lachnellula")+
  theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
         axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
         axis.title= element_text(face = "bold", size = 14, colour = "black"), 
         panel.background = element_blank(), legend.position="none",
         panel.border = element_rect(fill = NA, colour = "black")) + 
  ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) 

#################################################
# Kruskal Wallis to compare the abundance of different phyla between streams
###################################################

# agglomerate by genus
phylum.glom <- tax_glom(physeq.prop, taxrank = 'Phylum')

# get genus abundance matrix from phyloseq object
dat_phylum <- psmelt(phylum.glom)
dat_phylum$Phylum <- as.character(dat_phylum$Phylum)
dat_phylum$Phylum <- dat_phylum$Phylum %>% replace_na('Unclassified')
phylum.abund <- dat_phylum[,c(3,5,10)]

# choose phylum
phylum.final <- phylum.abund[which(phylum.abund$Phylum=="Basidiobolomycota"),]

# but really need to run Kruskas Wallis test for comparing more than two groups in order
# not to increase Type I error of multiple t-tests
kruskal.test(Abundance ~ Stream, data = phylum.final)

# then use pairwise comparisons between group levels with corrections for multiple testing
pairwise.wilcox.test(phylum.final$Abundance, phylum.final$Stream, p.adjust.method = "BH",
                     exact=F)

tgc <- summarySE(phylum.final, measurevar="Abundance", groupvars=c("Stream"))
ggplot(tgc, aes(x=Stream, y=Abundance, fill=Stream)) + geom_bar(stat="identity") + 
  theme_classic() + ggtitle("Mucoromycota")+
  theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
         axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
         axis.title= element_text(face = "bold", size = 14, colour = "black"), 
         panel.background = element_blank(), legend.position="none",
         panel.border = element_rect(fill = NA, colour = "black")) +
    ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,
                position=position_dodge(.9)) 

# this function writes p value of Kruskal Wallis test above chart, but to use this, need
# to use original data and not summarized
stat_compare_means()
ggbarplot(phylum.final, x="Stream", y="Abundance", add="mean_se") +stat_compare_means() 

geom_bar(position = "dodge", 
         stat = "summary", fun = "mean")

#################################################
# Kruskal Wallis test to compare the abundance of different phylum by transect at Hansen
# Happy, and Yako
###################################################

# choose stream
physeq_nobog <- subset_samples(physeq.prop, Transect != "Bog")
physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
physeq_happy <- subset_samples(physeq_nobog, Stream == "Happy")
physeq_yako <- subset_samples(physeq.prop, Stream == "Yako")

# agglomerate by genus
phylum.glom <- tax_glom(physeq_yako, taxrank = 'Phylum')

dat_phylum <- psmelt(phylum.glom)
dat_phylum$Phylum <- as.character(dat_phylum$Phylum)
dat_phylum$Phylum <- dat_phylum$Phylum %>% replace_na('Unclassified')
phylum.abund <- dat_phylum[,c(3,6,10)]

# create unique number for each genera
phylum.abund <- transform(phylum.abund,id=as.numeric(factor(Phylum)))

# create vector for each genera and p value
pvalue <- c(1:10)
w1 <- c(1:10)
w2 <- c(1:10)
w3 <- c(1:10)

for (i in 1:10){
  phylum.final <- phylum.abund[which(phylum.abund$id==i),]
  k <- kruskal.test(Abundance ~ Transect, data = phylum.final)
  pvalue[i] <- k$p.value
  w <- pairwise.wilcox.test(phylum.final$Abundance, phylum.final$Transect, 
                            p.adjust.method = "BH", exact=F)
  w1[i] <- w$p.value[1,1]
  w2[i] <- w$p.value[2,1]
  w3[i] <- w$p.value[2,2]
}

# create dataframe with phylum and p value combined
phylum.p <- data.frame(levels(factor(phylum.abund$Phylum)), pvalue, w1, w2, w3)
phylum.p <- phylum.p[order(pvalue),]
phylum.sig <- phylum.p[which(phylum.p$pvalue < 0.051),]
#colnames(phylum.sig) <- c("Phylum","P value","S1_S4","S1_S6","S4_S6")
#colnames(phylum.sig) <- c("Phylum","P value","H1_H2","H1_H3","H2_H3")
colnames(phylum.sig) <- c("Phylum","P value","Y1_Y2","Y1_Y3","Y2_Y3")
phylum.sig <- phylum.sig[,-2]
View(phylum.sig)

# choose phylum
phylum.final <- phylum.abund[which(phylum.abund$Phylum=="Basidiomycota"),]

# but really need to run Kruskas Wallis test for comparing more than two groups in order
# not to increase Type I error of multiple t-tests
kruskal.test(Abundance ~ Transect, data = phylum.final)

# then use pairwise comparisons between group levels with corrections for multiple testing
pairwise.wilcox.test(phylum.final$Abundance, phylum.final$Transect, p.adjust.method = "BH",
                     exact=F)

# whoops should have been using Dunn's test
library(FSA)
library(dunn.test)
dunnTest(Abundance ~ Transect, data=phylum.final, method="bonferroni")

tgc <- summarySE(phylum.final, measurevar="Abundance", groupvars=c("Transect"))
ggplot(tgc, aes(Transect, Abundance, fill=Transect)) + geom_bar(position = "dodge", 
  stat = "summary", fun = "mean") + theme_classic() + ggtitle("Basidiomycota")+
  theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
         axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
         axis.title= element_text(face = "bold", size = 14, colour = "black"), 
         panel.background = element_blank(), legend.position="none",
         panel.border = element_rect(fill = NA, colour = "black")) + 
  ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) 


#################################################
# Wilcoxon test to compare the abundance of different phylum by Hansen bank
###################################################

# subset by Hansen and agglomerate by phylum
physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
physeq_hansen_close <- subset_samples(physeq_hansen, Distance < 7)
phylum.glom <- tax_glom(physeq_hansen_close, taxrank = 'Phylum')


# get genus abundance matrix from phyloseq object
dat_phylum <- psmelt(phylum.glom)
dat_phylum$Phylum <- as.character(dat_phylum$Phylum)
dat_phylum$Phylum <- dat_phylum$Phylum %>% replace_na('Unclassified')
phylum.abund <- dat_phylum[,c(3,7,10)]

# create unique number for each genera
phylum.abund <- transform(phylum.abund,id=as.numeric(factor(Phylum)))

# create vector for each genera and p value
pvalue <- c(1:10)

for (i in 1:10){
  phylum.final <- phylum.abund[which(phylum.abund$id==i),]
  left <- phylum.final$Abundance[which(phylum.final$Side=="SL")]
  right <- phylum.final$Abundance[which(phylum.final$Side=="SR")]
  w <- wilcox.test(left, right, exact=F)
  pvalue[i] <- w$p.value
}

# create dataframe with genus and p value combined
phylum.p <- data.frame(levels(factor(phylum.abund$Phylum)), pvalue)
phylum.p <- phylum.p[order(pvalue),]
phylum.sig <- phylum.p[which(phylum.p$pvalue < 0.051),]
colnames(phylum.sig) <- c("Phylum","P value")
View(phylum.sig)

# do each phylum manually
# choose phylum
phylum.final <- phylum.abund[which(phylum.abund$Phylum=="Rozellomycota"),]
left <- phylum.final$Abundance[which(phylum.final$Side=="SL")]
right <- phylum.final$Abundance[which(phylum.final$Side=="SR")]

# then use pairwise comparisons between group levels with corrections for multiple testing
wilcox.test(left, right, exact=F)

tgc <- summarySE(phylum.final, measurevar="Abundance", groupvars=c("Side"))
ggplot(tgc, aes(Side, Abundance, fill=Side)) + geom_bar(position = "dodge", 
       stat = "summary", fun = "mean") + theme_classic() + ggtitle("Rozellomycota")+
  theme( axis.text.x = element_text(face = "bold",colour = "black", size = 12), 
         axis.text.y = element_text(face = "bold", size = 11, colour = "black"), 
         axis.title= element_text(face = "bold", size = 14, colour = "black"), 
         panel.background = element_blank(), legend.position="none",
         panel.border = element_rect(fill = NA, colour = "black")) + 
  ylab("Relative Abundance") + 
  geom_errorbar(aes(ymin=Abundance-se, ymax=Abundance+se), width=.2,position=position_dodge(.9)) 
